PROGRAM TorsionTest ! Torsion of a square cross-section bar. ! Example taken from Section 8.4.3 on Solid Mechanics in the following: ! Reddy, J. N.: "An Introduction to the Finite Element Method", 2nd ed., McGraw-Hill . ! Copyright JGX Software Solutions, 2018. ! ! THIS SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN ! THIS SOFTWARE. ! Note that this source code is designed for use with the Intel(R) Math Kernel Library ! version of Pardiso. USE, INTRINSIC :: ISO_FORTRAN_ENV IMPLICIT NONE INTEGER, PARAMETER :: dp = REAL64 INTEGER(INT64), DIMENSION(64) :: pt ! Pardiso input data. INTEGER, DIMENSION(64) :: iparm ! See https://software.intel.com/en-us/mkl-developer-reference-fortran-pardiso-iparm-parameter INTEGER mtype ! and https://software.intel.com/en-us/mkl-developer-reference-fortran-pardiso INTEGER maxfct, mnum, phase, n ! for details. INTEGER nrhs, msglvl ! Pardiso output data. INTEGER error INTEGER perm(64) ! Space for the permutation vector calculated by Pardiso. REAL(REAL64), ALLOCATABLE :: a(:) ! Space to hole the matrix of FD coefficients as well as INTEGER, ALLOCATABLE :: ia(:), ja(:) ! the matrix's topology in 3-array CRS format. REAL(REAL64), ALLOCATABLE :: b(:,:) ! Space for the RHS. REAL(REAL64), ALLOCATABLE :: x(:,:) ! Space for the solution. REAL(REAL64) psi(0:7,0:7) ! An explicit interface is required for InitializeTorsion8x8 as it ! has a dummy argument with the ALLOCATE attribute. INTERFACE SUBROUTINE InitializeTorsion8x8(a, ia, ja) USE, INTRINSIC :: ISO_FORTRAN_ENV REAL(REAL64), ALLOCATABLE :: a(:) INTEGER, ALLOCATABLE :: ia(:), ja(:) END SUBROUTINE InitializeTorsion8x8 END INTERFACE !------------------------------------------------------------------------------ ALLOCATE(x(64,1), b(64,1)) mtype = 11 ! defines a real and non-symmetric matrix. ! Pardiso initialization. CALL pardisoinit(pt, mtype, iparm) maxfct = 1 ! We are working with just one matrix. mnum = 1 ! We are factorizing the first of just one matrix ! phase = 13 ! We require analysis, numerical factorization, solve, iterative refinement. n = 64 ! Solving for 64 equations in Ax = b. ! Grab the matrix a in 3-array CSR format. CALL InitializeTorsion8x8(a, ia, ja) nrhs = 1 ! Solving for just one right hand side. msglvl = 0 ! Set the message level to none. ! Initialize the RHS. b = -2.0_dp / 256.0_dp iparm(5) = 2 ! Return the permutation vector computed during phase 1. iparm(60) = 0 ! Use in-core (IC) mode. If memory problems occur, set to 2, out-of-core (OOC) mode. iparm(27) = 0 ! Matrix checker set to 0. Set to 1 in order to diagnose problems. ! Solve for the stress function values. On exit from the following routine, the stress function ! values will be placed in X. CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, perm, nrhs, iparm, msglvl, b, x, error) ! The stress function values can also be indexed by the original (i,j) values .. psi = TRANSPOSE(RESHAPE(X(1:64,1), (/ 8, 8 /))) END PROGRAM TorsionTest SUBROUTINE InitializeTorsion8x8(a, ia, ja) ! Copyright JGX Software Solutions, 2018. ! ! THIS SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN ! THIS SOFTWARE. ! Formulate a sparse 64x64 matrix of finite difference coefficients in 3-array CSR format. ! The matrix is taken from a finite difference approximation for the shear stresses incurred ! due to the torision of a solid square bar of side length 1 unit. ! See Example 8.8 from J.N.Reddy, "Finite Element Method" for details. ! For details on 3-array CSR format see ! https://software.intel.com/en-us/mkl-developer-reference-fortran-sparse-blas-csr-matrix-storage-format#CSR3 USE, INTRINSIC :: ISO_FORTRAN_ENV IMPLICIT NONE REAL(REAL64), ALLOCATABLE :: a(:) INTEGER, ALLOCATABLE :: ia(:), ja(:) INTEGER, PARAMETER :: dp = REAL64 INTEGER ioff INTEGER K !------------------------------------------------------------------------------ ALLOCATE(a(288)) ALLOCATE(ia(65), ja(288)) ! IA(ir) contains the element of JA that contains the first column index of row ir. ia(1:8) = (/ 1, 4, 8, 12, 16, 20, 24, 28 /) ja(1:7) = (/ 4, 9, 14, 19, 24, 29, 34 /) ! ja used as temporary storage. ioff = 31 DO K = 1, 6 ia(K*8+1) = ioff ia(K*8+2 : (K+1)*8) = ioff + ja(1:7) ioff = ia((K+1)*8) + 4 END DO ioff = ia(56) + 3 ia(57:64) = ia(1:8) + ioff ia(65) = ia(64) + 3 a(1:30) = (/ -4.0_dp, 2.0_dp, 2.0_dp, & 1.0_dp, -4.0_dp, 1.0_dp, 2.0_dp, & 1.0_dp, -4.0_dp, 1.0_dp, 2.0_dp, & 1.0_dp, -4.0_dp, 1.0_dp, 2.0_dp, & 1.0_dp, -4.0_dp, 1.0_dp, 2.0_dp, & 1.0_dp, -4.0_dp, 1.0_dp, 2.0_dp, & 1.0_dp, -4.0_dp, 1.0_dp, 2.0_dp, & 1.0_dp, -4.0_dp, 2.0_dp /) a(31:68) = (/ 1.0_dp, -4.0_dp, 2.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp /) a(69:106) = a(31:68) a(107:144) = a(31:68) a(145:182) = a(31:68) a(183:220) = a(31:68) a(221:258) = a(31:68) a(259:288) = (/ 1.0_dp, -4.0_dp, 2.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp, 1.0_dp, & 1.0_dp, 1.0_dp, -4.0_dp /) ! JA(iv) contains the column index of value iv in a. ja(1:30) = (/ 1, 2, 9, & 1, 2, 3, 10, & 2, 3, 4, 11, & 3, 4, 5, 12, & 4, 5, 6, 13, & 5, 6, 7, 14, & 6, 7, 8, 15, & 7, 8, 16 /) ja(31:68) = (/ 1, 9, 10, 17, & 2, 9, 10, 11, 18, & 3, 10, 11, 12, 19, & 4, 11, 12, 13, 20, & 5, 12, 13, 14, 21, & 6, 13, 14, 15, 22, & 7, 14, 15, 16, 23, & 8, 15, 16, 24 /) ja(69:106) = ja(31:68) + 8 ja(107:144) = ja(31:68) + 16 ja(145:182) = ja(31:68) + 24 ja(183:220) = ja(31:68) + 32 ja(221:258) = ja(31:68) + 40 ja(259:288) = (/ 1, 9, 10, & 2, 9, 10, 11, & 3, 10, 11, 12, & 4, 11, 12, 13, & 5, 12, 13, 14, & 6, 13, 14, 15, & 7, 14, 15, 16, & 8, 15, 16 /) ja(259:288) = ja(259:288) + 48 END SUBROUTINE InitializeTorsion8x8